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Abstract 

The paper presents a general theory of couphng of eigenvalues of complex matrices of 
arbitrary dimension depending on real parameters. The cases of weak and strong coupling are 
distinguished and their geometric interpretation in two and three-dimensional spaces is given. 
General asymptotic formulae for eigenvalue surfaces near diabolic and exceptional points are 
presented demonstrating crossing and avoided crossing scenarios. Two physical examples 
illustrate effectiveness and accuracy of the presented theory. 

PACS numbers: 02.10.Yn, 02.30.Oz, 42.25.Bs 

1 Introduction 

Behavior of eigenvalues of matrices dependent on parameters is a problem of general interest hav- 
ing many important applications in natural and engineering sciences. Probably, [Hamilton (1833)| 



was the first who revealed an interesting physical effect associated with coincident eigenvalues 
known as conical refraction, see also 



Berry et al. (1999) 



In modern physics, e.g. quantum 
mechanics, crystal optics, physical chemistry, acoustics and mechanics, singular points of ma- 
trix spectra associated with specific effects attract great interest of researchers since the papers 
[Von Neumann and Wigner (1929) Herring (1937)[ Teller (1937)| . These are the points where 



matrices possess multiple eigenvalues. In applications the case of double eigenvalues is the most 
important. With a change of parameters coupling and decoupling of eigenvalues with crossing 
and avoided crossing scenario occur. The crossing of eigenvalue surfaces (energy levels) is con- 



nected with the topic of geometrical phase, see Berry and Wilkinson (1984) 


. In recent papers, see 


e.g. [Heiss (2000) 


Dembowsky et al. (2001) Berry and Dennis (2003)| Dembowsky et al. (2003)| 


Keck et al. (2003) 


Korsch and Mossman (2003)| Heiss (2004)| 


Stehmann et al. (2004)|, two im- 



portant cases are distinguished: the diabolic points (DPs) and the exceptional points (EPs). From 
mathematical point of view DP is a point where the eigenvalues coalesce (become double), while 
corresponding eigenvectors remain different (linearly independent); and EP is a point where both 
eigenvalues and eigenvectors merge forming a Jordan block. Both the DP and EP cases are in- 
teresting in applications and were observed in experiments, see e.g. Dembowsky et al. (2001)| 



Dembowsky et al. (2003) Stehmann et al. (2004) . In early studies only real and Hermitian ma- 
trices were considered while modern physical systems require study of complex symmetric and non- 
symmetric matrices, see Mondragon and Hernandez (1993) Berry and Dennis (2003) Keck et al. (2003) 



Note that most of the cited papers dealt with specific 2x2 matrices depending on two or three 
parameters. Of course, in the vicinity of an EP (and also DP) the m-dimensional matrix problem 
becomes effectively two-dimensional, but finding the corresponding two-dimensional space for a 
general m-dimensional matrix family is a nontrivial problem [Arnold (1983)| . 
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In this paper we present a general theory of couphng of eigenvalues of complex matrices of arbi- 
trary dimension smoothly depending on multiple real parameters. Two essential cases of weak and 
strong coupling based on a Jordan form of the system matrix are distinguished. These two cases 
correspond to diabolic and exceptional points, respectively. We derive general formulae describing 
coupling and decoupling of eigenvalues, crossing and avoided crossing of eigenvalue surfaces. We 
present typical (generic) pictures showing movement of eigenvalues, the eigenvalue surfaces and 
their cross-sections. It is emphasized that the presented theory of coupling of eigenvalues of com- 
plex matrices gives not only qualitative, but also quantitative results on behavior of eigenvalues 
based only on the information taken at the singular points. Two examples on propagation of light 
in a homogeneous non-magnetic crystal possessing natural optical activity (chirality) and dichro- 
ism (absorption) in addition to biaxial birefringence illustrate basic ideas and effectiveness of the 
developed theory. 

The presented theory is based on previous research on interaction of eigenvalues of real ma- 
trices depending on multiple parameters with mechanical applications. In Seyranian (1991)] 



Seyranian (1993) the important notion of weak and strong coupling (interaction) was intro- 
duced for the first time. In the papers [Seyranian and Pedersen (1993)[ Seyranian et al. (1994) 



Mailybaev and Seyranian (1999) Seyranian and Kliem (2001)| Seyranian and Mailybaev (2001) 



Kirillov and Seyranian (2002)||Seyranian and Mailybaev (2Q03)||Kirillov (2004)[|Kirillov and Seyranian (2004) 
and in the recent book [Seyranian and Mailybaev (2003)| significant mechanical effects related to 
diabolic and exceptional points were studied. These include transference of instability between 
eigenvalue branches, bimodal solutions in optimal structures under stability constraints, flutter 
and divergence instabilities in undamped nonconservative systems, effect of gyroscopic stabiliza- 
tion, destabilization of a nonconservative system by infinitely small damping, which were described 
and explained from the point of view of coupling of eigenvalues. An interesting application of the 
results on eigenvalue coupling to electrical engineering problems is given in Dobson et al. (2001)[ . 



The paper is organized as follows. In Section 2 we present general results on weak and strong 
coupling of eigenvalues of complex matrices depending on parameters. These two cases correspond 
to the study of eigenvalue behavior near diabolic and exceptional points. Section 3 is devoted 
to crossing and avoided crossing of eigenvalue surfaces near double eigenvalues with one and two 
eigenvectors. Two physical examples are presented in Section 4, and finally we end up with the 
conclusion in Section 5. 



2 Coupling of eigenvalues 

Let us consider the eigenvalue problem 

Au = Au (1) 

for a general m x m complex matrix A smoothly depending on a vector of n real parameters p = 
(pi, . . . ,Pn)- Assume that, at p = po, the eigenvalue coupling occurs, i.e., the matrix Aq ~ A(po) 
has an eigenvalue Aq of multiplicity 2 as a root of the characteristic equation det(Ao — AqI) = 0; I is 
the identity matrix. This double eigenvalue can have one or two linearly independent eigenvectors 
u, which determine the geometric multiplicity. The eigenvalue problem adjoint to is 

A*v = 77V, (2) 



where A* — A^ is the adjoint matrix operator (Hermitian transpose), see e.g. Lancaster (1969)| . 



The eigenvalues A and 77 of problems (Q) and J^J are complex conjugate: 77 = A. 

Double eigenvalues appear at sets in parameter space, whose codimensions depend on the matrix 
type and the degeneracy (EP or DP). In Tabled we list these codimensions based on the results of 



the singularity theory Von Neumann and Wigncr (1929)j,Arnold (1983), . In this paper we analyze 



general (nonsymmetric) complex matrices. The EP degeneracy is the most typical for this type of 
matrices. In comparison with EP, the DP degeneracy is a rare phenomenon in systems described 
by general complex matrices. However, some nongeneric situations may be interesting from the 
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matrix type 


codimcnsion of DP 


codimcnsion of EP 


real symmetric 


2 


non-existent 


real nonsymmctric 


3 


1 


Hcrmitian 


3 


non-existent 


complex symmetric 


4 


2 


complex nonsymmctric 


6 


2 



Table 1: Codiniensions of eigenvalue degeneracies. 



physical point of view. As an example, we mention complex non-Hermitian perturbations of 
symmetric two-parameter real matrices, when the eigenvalue surfaces have coffee-filter singularity, 
Mondragon and Hernandez (1993){ Berry and Dennis (2003)} Keck et al. (2003)]. A general 



see 



theory of this phenomenon will be given in our companion paper [Kirillov et al. (2004)| 

Let us consider a smooth perturbation of parameters in the form p = p(e), where p(0) = Po 
and e is a small real number. For the perturbed matrix A = A(p(e)), we have 



A = Ao + eAi + ie^Aa + o{e^), 

A Af ) A A '^^^^ 'V^ ^^"^ 

Po , 1 2^ Q , 2 2^ Q 2^ dpidpi de de 



(3) 



The double eigenvalue Aq generally splits into a pair of simple eigenvalues under the perturba- 
tion. Asymptotic formulae for these eigenvalues and corresponding eigenvectors contain integer or 
fractional powers of e [Vishik and Lyusternik (1960) . 



2.1 Weak coupling of eigenvalues 

Let us consider the coupling of eigenvalues in the case of Ao with two linearly independent eigen- 
vectors Ui and U2. This coupling point is known as a diabolic point. Let us denote by vi and V2 
two eigenvectors of the complex conjugate eigenvalue rj = X for the adjoint eigenvalue problem ^ 
satisfying the normalization conditions 

(Ui, Vi) = (U2, V2) 1, (Ui, V2) = (U2, Vi) = 0, (4) 

where (u, v) = X]"=i '^i^i denotes the Hermitian inner product. Conditions Q define the unique 
vectors Vi and V2 for given Ui and U2 Seyranian and Mailybaev (2003)| . 



For nonzero small e, the two eigenvalues A+ and A_ resulting from the bifurcation of Ao and 
the corresponding eigenvectors u± are given by 

X± ^ Xo + ^i±e + o{e), u± = a±ui -f /3±U2 -I- o(l). (5) 

The coefficients fj,±, a±, and P± are found from the 2x2 eigenvalue problem 
(see e.g. Seyranian and Mailybaev (2003)| ) 



(AiUi,vi) (AiU2,vi) \ / a± \ _ / a± 

(Aiui,v2) (AiU2,V2) \ P± J ^^^\ I3± ^ ' 



Solving the characteristic equation for we find 



(AiUi,Vi) + (AiU2,V2) /((AiUi,Vi) - (AiU2,V2))2 

f^± = ^ ±Y ^ -f (AiUi,V2)(AiU2,vi). (7) 

We note that for Hermitian matrices A one can take Vi = Ui and V2 = U2 in ©, where the 
eigenvectors Ui and U2 are chosen satisfying the conditions (ui, Ui) = (u2, U2) = 1 and (ui, U2) = 
0, and obtain the well-known formula, see [Courant and Hilbert (1953)] . 
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Figure 1: Eigenvalue coupling: (a) weak, (b) strong. 



As the parameter vector passes the coupling point po along the curve p(£) in parameter space, 
the eigenvalues A+ and A„ change smoothly and cross each other at Aq, see Figure^. At the 
same time, the corresponding eigenvectors u+ and u_ remain different (linearly independent) at all 
values of s including the point po. We call this interaction weak coupling. By means of eigenvectors, 
the eigenvalues A± are well distinguished during the weak coupling. 

We emphasize that despite the eigenvalues \± and the eigenvectors u± depend smoothly on a 
single parameter e, they are nondifferentiable functions of multiple parameters at po in the sense 



ofFrechet Schwartz (1967) 



2.2 Strong coupling of eigenvalues 

Let us consider coupling of eigenvalues at po with a double eigenvalue Ao possessing a single eigen- 
vector Uq. This case corresponds to the exceptional point. The second vector of the invariant 
subspace corresponding to Ao is called an associated vector Ui (also called a generalized eigenvec- 
tor [Lancaster (1969) ); it is determined by the equation 



AqUi AqUi -hUQ. (8) 

An eigenvector vq and an associated vector vi of the matrix A* are determined by 

AqVq Aqvo, AqVi = Aqvi + vo, (ui,vo)==1, (ui,vi)==0, (9) 

where the last two equations are the normalization conditions determining vo and vi uniquely for 
a given Ui. 

Bifurcation of Aq into two eigenvalues A± and the corresponding eigenvectors u± are described 



by (see e.g. Seyranian and Mailybaev (2003) ) 



A± Ao ± ^/mTe + Ai2e + o(e), ^^^^ 
u± = Uo ± Ui^/Til? + ((UiUq -I- ^2Ui - G~^AiUo)e + o(e), 

where G = Aq — Aol + UivJ;. The coefficients /ii and ^2 are 

^1 = (AiUo, Vo), ^2 ((AiUo, vi) -f (AiUi, Vo))/2. (11) 

With a change of e from negative to positive values, the two eigenvalues A± approach, collide 
with infinite speed (derivative with respect to e tends to infinity) at Ao, and diverge in the per- 
pendicular direction, see Figure QJj. The eigenvectors interact too. At e = 0, they merge to Uo 
up to a scalar complex factor. At nonzero £, the eigenvectors u± differ from Uo by the leading 
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term ±Ui^/iie. This term takes the purely imaginary factor i as e changes the sign, for example 
altering from negative to positive values. 

We call such a coupling of eigenvalues as strong. An exciting feature of the strong coupling is 
that the two eigenvalues cannot be distinguished after the interaction. Indeed, there is no natural 
rule telling how the eigenvalues before coupling correspond to those after the coupling. 

3 Crossing of eigenvalue surfaces 

3.1 Double eigenvalue with single eigenvector 

Let, at the point po, the spectrum of the complex matrix family A(p) contain a double complex 
eigenvalue Aq with an eigenvector Uq and an associated vector Ui. The splitting of the double 
eigenvalue with a change of the parameters is governed by equations (fTn|l and (fTT|l . Introducing 
the real n-dimensional vectors f , g, h, r with the components 




s = 1, . . . , n. 



and neglecting higher order terms, we obtain from (|10|l an asymptotic formula 



ReAA + ilmAA = ±^{i, Ap> + i(g, Ap> + -((h, Ap) + i(r, Ap)), (14) 

where A A Aj- — Ao, Ap = p — Po, and angular brackets denote inner product of real vectors: 
(a, b) = X]l=i '^sbs- From equation (|14|l it is clear that the eigenvalue remains double in the first 
approximation if the two following equations are satisfied 

(f,Ap)=0, (g,Ap)=0. (15) 

This means that the double complex eigenvalue with the Jordan chain of length 2 has codimension 
2. Thus, double complex eigenvalues occur at isolated points of the plane of two parameters, and 
in the three-parameter space the double eigenvalues form a curve [Arnold (1983), . Equations H15|l 
define a tangent line to this curve at the point po- 

Taking square of (|14|l . where the terms linear with respect to the increment of parameters are 
neglected, and separating real and imaginary parts, we derive the equations 

(ReAA)2 - (ImAA)2 = (f, Ap), 2RcAAImAA = (g, Ap). (16) 

Let us assume that fi+gfy^ 0, which is the nondegeneracy condition for the complex eigenvalue 
Aq. Isolating the increment Api in one of the equations H16|l and substituting it into the other one 
we get 

gi{ReAXf - 2/iReAAImAA - gi{lmAX)^ = 7, (17) 
where 7 is a small real constant 

n 

7^J2^fsgi-figs)Ap,. (18) 

s=2 

Equation (|17|l describes hyperbolic trajectories of the eigenvalues X± in the complex plane when 
only Api is changed and the increments Ap2, . . ., Ap„ are fixed. Of course, any component of the 
vector Ap can be chosen instead of Api . 
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Let us study movement of eigenvalues in the complex plane in more detail. If Apj — 0, 
j = 2, . . . , n, or if they are nonzero but satisfy the equality 7 = 0, then equation H17|l yields two 
perpendicular lines which for gi ^ are described by the expression 

5iRe(A - Ao) - (^/i ± Ini(A - Aq) = 0. (19) 

These lines intersect at the point Aq of the complex plane. Due to variation of the parameter pi 
two eigenvalues A± approach along one of the lines (|19|l . merge to Aq at Api = 0, and then diverge 
along the other line l|19|l . perpendicular to the line of approach; see Figure Eb, where the arrows 
show motion of eigenvalues with a monotonous change of pi . Recall that the eigenvalues that born 
after the coupling cannot be identified with the eigenvalues before the coupling. 

If 7 0, then equation (|17|l defines a hyperbola in the complex plane. Indeed, for gi 7^ it is 
transformed to the equation of hyperbola 

(giRe(A - Aq) - /iIm(A - Ao))^ - (Im(A ~ X„)f{f? + gl) = 751 (20) 

with the asymptotes described by equation (|19|l . As Api changes monotonously, two eigenvalues 
A+ and A_ moving each along its own branch of hyperbola come closer, turn and diverge; see 
Figure |3i,c. Note that for a small 7 the eigenvalues A± come arbitrarily close to each other 
without coupling that means avoided crossing. When 7 changes the sign, the quadrants containing 
hyperbola branches are changed to the adjacent. 

Expressing ImAA from the second of equations H16|l . substituting it into the first equation and 
then isolating ReAA, we find 



ReA± = Ao + -(h, Ap) ± ^- ((f , Ap) + ^(f , Ap)^ + (g, Ap)^) . (21) 
Similar transformation yields 



ImA± = Ao + - (r , Ap) ± ^ - (- (f , Ap) + ^(f , Ap)^ + (g, Ap)^) . (22) 

Equations (|21|l and (|22|l describe behavior of real and imaginary parts of eigenvalues A± 
with a change of the parameters. On the other hand they define hypersurfaces in the spaces 
(pijP2, • ■ • jPn,ReA) and (pi,?32, ■ • ■ jPmlmA). The sheets ReA+(p) and ReA_(p) of the eigenvalue 
hypersurface H21|l are connected at the points of the set 

ReAA ^ i (h, Ap) , (g, Ap) = 0, (f , Ap) < 0, (23) 

where the real parts of the eigenvalues A± coincide: ReA_ = RcA+. Similarly, the set 

ImAA = i(r, Ap), (g, Ap) = 0, (f , Ap) > 0, (24) 

glues the sheets ImA+(p) and ImA_(p) of the eigenvalue hypersurface (|22|l . 

To study the geometry of the eigenvalue hypersurfaces we look at their two-dimensional cross- 
sections. Consider for example the functions ReA(pi) and ImA(pi) at fixed values of the other 
parameters p2,P3, ■ ■ • ,Pn- When the increments Aps = 0, s = 2,3, ...,n, both the real and 
imaginary parts of the eigenvalues A± cross at pi = Pi, see Figure The crossings are described 
by the double cusps defined by the equations following from 1)21(1 and ((22|l as 

ReAA.±^A±A^E^I^+^Ap„ ImAA = ±^^A±^Z?±ZA,, + ^A,, (25) 
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Figure 3: Crossing of eigenvalue surfaces near the double eigenvalue with single eigenvector. 



For the fixed Aps 7^ 0, s — 2, 3, . . . ,n, either real parts of the eigenvalues A± cross due to 
variation of pi while the imaginary parts avoid crossing or vice-versa, as shown in Figure [2^,c. 
Note that these two cases correspond to level crossing and width repulsion or vice- versa studied in 
Heiss (2000)1 ■ The crossings, which occur at = Pi ~ ^^^2^9^/ Qij^Vs a^nd 



^ n 1 

ReA/i = ReAo - — V'(ft.i.gs - gihs)Aps, ImA^ = ImAo - — y^(rigs - girs)Aps, (26) 

are described by the equations ()21|l and H22|l . In the vicinity of the crossing points the tangents of 
two intersecting curves are 

ReA = ReA. + (^^ ± [p, -p^), (27) 



ImA = ImA, + (^!^ ± j (p, - ), (28) 

where the coefficient 7 is defined by equation 118|l . Lines H27|) and H28|l tend to the vertical position 
as 7 — > and coincide at 7 = 0. The avoided crossings are governed by the equations (|21|l and 

(EH). 

If the vector of parameters consists of only two components p = (pi,p2), then in the vicinity of 
the point po, corresponding to the double eigenvalue Aq, the eigenvalue surfaces H21|) and H22|) have 
the form of the well-known Whitney umbrella; see Figure |21 The sheets of the eigensurfaces are 
connected along the rays H23|l and (|24|l . We emphasize that these rays are inclined with respect to 
the plane of the parameters pi, p2- The cross-sections of the eigensurfaces by the planes orthogonal 
to the axis p2, described by the equations are shown in Figure 12 

Note that the rays (|23|l . (|24|l and the point po are well-known in crystal optics as the branch 
cuts and the singular axis, respectively Berry and Dennis (2003)] . We emphasize that the branch 



cut is a general phenomenon: it always appears near the EP degeneracy. In general, branch cuts 
may be infinite or end up at another EP. The second scenario is always the case when the complex 
matrix A(p) is a small perturbation of a family of real symmetric matrices Kirillov et al. (2004)| . 



Consider the movement of the eigenvalues in the complex plane near the point po due to cyclic 
variation of the parameters pi and p2 of the form Api = a + r cos tp and Ap2 = b + r sin ip, where 
a, 6, and r are small parameters of the same order. From equations H16|) we derive 

{giReAX" - 2/iReAAImAA - gilmAA^ - b{f2gi - hg2)f+ 
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a) 



Im X 



Im X 



b) 




Hex Q 



Im X 



c) 



d) 




Re 



Im X 




ReX 



Im X 



Re X 



Figure 4: Movement of eigenvalues due to cyclic evolution of the parameters. 



(.gaRcAA^ - 2/2ReAAImAA - galmAA^ - a{hg2 - 51/2))' = (/2gi ~ h92fr^ 



(29) 



Movement of eigenvalues on the complex plane governed by equation (|29|l is shown in Figure^ 
If the contour encircles the point po, then the eigenvalues move along the curve H29|l around the 
double eigenvalue Aq in the complex plane, see Figure^. Indeed, in this case a? A-b"^ < and the 
loop (|29|l crosses the lines ReA = ReAo and ImA = ImAo at the four points given by the equations 



(ImAA)2 = 



(/2.gi - /iff2) (ff2a - gib ± \/ ima - gibf + {t 



9l 



9l 



and 



(ReAA)2 



(/251 " fm) {gib - 52a ± Vi9ib - g2ay + (r^ ~ - b^){gj + 



5? + 5i 



(30) 



(31) 



respectively. When + = the loop overlaps at the double eigenvalue and its form depends 
on the sign of the quantity a — (/251 — fig2){gib — 320). If cr < the eigenvalues cross the line 
ReA — ReAo (Figure ^Is), otherwise they cross the line ImA — ImAo (Figure 0Ji). Eigenvalues 
strongly couple at the point Aq in the complex plane. For + b^ > the circuit in the parameter 
plane does not contain the point po and the eigenvalues move along the two different closed paths 
("kidneys" , [Arnold (1989)| ) in the complex plane, see Figure0Ji,e. Each eigenvalue crosses the line 
ReA — ReAo twice for tr < (Figure0^), and for ct > they cross the line ImA — ImAo (Figure^). 
Note that the "kidneys" in the complex plane were observed in Korsch and Mossman (2003)| for 
the specific problem of Stark resonances for a double 6 quantum well. 

3.2 Double eigenvalue with two eigenvectors 

Let Ao be a double eigenvalue of the matrix Ao = A(po) with two eigenvectors Ui and U2. Under 
perturbation of parameters p = Po + Ap, the bifurcation of Ao into two simple eigenvalues A+ 
and A_ occurs. Using and Q, we obtain the asymptotic formula for X± under multiparameter 
perturbation as 



Ao 



(dii 



.Ap) 



± 



;dii -d22,Ap)2 



(di2,Ap>(d21,Ap), 



where dy — (dj, , . . . , c?"- ) is a complex vector with the components 



. ^ (dA 
\dpk 



Ui, V, 



(32) 



(33) 
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and (dy, Ap) — (Redy, Ap) + z(Imdij, Ap). In the same way as we derived formulae (|21|) and 
(|22|l . we obtain from H32|l the expressions for real and imaginary parts of \± in the form 



ReA± =ReAo + Re(dn +d22,Ap)/2± v/(|c| +Rec)/2, (34) 



Im A± = Im Ao + Im (dn + daa, Ap)/2 ± y/{\c\ -Rec)/2, (35) 

where 

c= (dn-d22,Ap)V4+(di2,Ap)(d2i,Ap). (36) 

Considering the situation when Ao remains double under perturbation of parameters, i.e. A+ = 
A_, we obtain the two independent equations 

Rec = 0, lmc = 0. (37) 

By using ©-ijTJ), one can show that the perturbed double eigenvalue A4. — A_ possesses a 
single eigenvector u_|_ — u_, i.e., the weak coupling becomes strong due to perturbation, see 
[Seyranian and Mailybaev (2003) 



The perturbed double eigenvalue has two eigenvectors only when the matrix in the left-hand 
side of ^ is proportional to the identity matrix. This yields the equations 

(dn , Ap) = (d22 , Ap) , (di2 , Ap) (d2i , Ap) ^ 0, (38) 

Conditions H38|) imply H37I) and represent six independent equations taken for real and imaginary 
parts. Thus, weak coupling of eigenvalues is a phenomenon of codimension 6, which generically oc- 



curs at isolated points in 6-parameter space, see [Arnold (1983) Mondragon and Hernandez (1993) 



First, let us study behavior of the eigenvalues A+ and A_ depending on one parameter, say 
pi , when the other parameters p2, ■ ■ ■ ,Pn are fixed in the neighborhood of the coupling point 
•^+(Po) = -^-(Po) = Ao. In case Ap2 = • • • = Ap„ = 0, expression (p?^ yields 



A± = Ao + i^Al±^ ± y + d}24i j Api. (39) 

The two eigenvalues couple when Api = with the double eigenvalue Ao, see Figure Et- As we 
showed in Section 2, the eigenvalues A-f and A_ behave as smooth functions at the coupling point; 
they possess different eigenvectors, which are smooth functions of Api too. 

If the perturbations Ap2, ■ • ■ , Ap„ are nonzero, the avoided crossing of the eigenvalues \± with 
a change of pi is a typical scenario. We can distinguish different cases by checking intersections of 
real and imaginary parts of A_)_ and A_. By using H34|l . we find that ReA_|_ = Re A_ if 

lmc = 0, Rec<0. (40) 

Analogously, from H35|l it follows that Im A+ = Im A- if 

lmc = 0, ReoO. (41) 

Let us write expression ()35f) in the form 

c = co + ciApi +c2(Api)2, (42) 



where 



Co 



= E [(4i-'^22)(4i-42)/4 + 424i] ApfcAp,, 



k,l=2 



ci = E [('^ii - 42)(4i - 42)/2 + idl24i + dUhl)] Apfc, 



(43) 



fc=2 



C2^idi,~di,)y4 + di,di,. 
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Figure 5: Weak coupling of eigenvalues and avoided crossing. 
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If the discriminant D = (Imci)^ — 4ImcoImc2 > 0, the equation Imc — yields two solutions 

^1 2IniC2 ^ 2Imc2 ^ ^ 

There are no real solutions if D < 0, and the single solution corresponds to the degenerate case 
D = 0. At the points Pi—Pi + Ap" and p\ = + ApJ the values of c are real, and we denote 
them by Ca and Cf,, respectively. According to (|40|l and (|41|l . the sign of Ca,b determines whether 
the real or imaginary parts of A± coincide at p^'^. 

In the nondegenerate case D ^ 0, there are four types of avoided crossing shown in Figure|Sl3~e. 
The first case corresponds to Z? < when both real and imaginary parts of the eigenvalues A± 
are separate at all pi, see Figure Eb- In other cases D > 0, so that there are two separate points 
Pi and p\. For the second type we have Ca > and Cb < 0, when both real and imaginary parts 
of X± have a single intersection, see Figure The equivalent situation when Ca < and Cf, > 
is obtained by interchanging the points pi and p\ in Figure [S];. The third type is represented 
by Ca.b < 0, when the real parts of X± have two intersections and ImA± do not intersect, see 
Figure [Sji- Finally, if Ca^b > 0, when the real parts of X± do not intersect and Im X± intersect 
at both p1 and p\ , see Figure The last column in Figure O shows behavior of the eigenvalues 
A± on the complex plane. In each of the cases b-e, the trajectories of eigenvalues on the complex 
plane may intersect and/or self- intersect, which can be studied by using expression H32|l . Note 
that intersections of the eigenvalue trajectories on the complex plane do not imply eigenvalue 
coupling since the eigenvalues A+ and A_ pass the intersection point at different values of pi . The 
small loops of the eigenvalue trajectories on the complex plane, shown in Figure|33,e, shrink as the 
perturbations of the parameters Ap2, Aps, . . . , Ap„ tend to zero. Finally, we mention that the case 
of Figure is the only avoided crossing scenario when the eigenvalues follow the initial directions 
on the complex plane after interaction. In the other three cases (b,d, and e) the eigenvalues 
interchange their directions due the interaction. 

Let us consider a system depending on two parameters pi and p2 with the weak coupling of 
eigenvalues at pi = pi and P2 = P2- The double eigenvalue Aq bifurcates into a pair X± under 
perturbation of the parameters Api and Ap2- Conditions (|40|l and H41|) determine the values of 
parameters, at which the real and imaginary parts of A± coincide. 

Let us write expression (I36|) in the form 

C = Cii{Apif + Ci2ApiAp2 + C22(Ap2)^ (45) 

where 

cn = Ki-dL)V4 + d}2rf2i, C22^idj,-42?/4 + dl24„ 

ci2 - Ki - dDKi - dh)/^ + d'udl, + dl^d\^. 

If the discriminant D' = (Imci2)^ — 41mciilmc22 > 0, the equation Imc = yields the two 
crossing lines 

la-. 2ImciiApi + (Imci2 + Ap2 = 0, 



(47) 

lb-. 2ImciiApi + (Imci2 - V-D')Ap2 = 0. 

There are no real solutions if D' < 0, and the lines la and 4 coincide in the degenerate case D' = 0. 
On the lines la,b the values of c are real numbers of the same sign; we denote 7a = signc for the 
line la, and 7;, = signc for the line lb- According to H4Q(I and (|41|l . the real or imaginary parts of 
A± coincide at laM for negative or positive 7a, b, respectively. 

One can distinguish four types of the graphs for Re A-|-(pi,p2) and ImA±(pi,p2) shown in 
Figure In nondegenerate case D' 7^ 0, the eigenvalues A+ and A_ are different for all parameter 
values except the initial point pi^2 = Pi 2- If ^' ^ ^^e eigenvalue surfaces are cones with non- 
elliptic cross-section, see Figure [SK- Other three types correspond to the case D' > 0. If 7a < 
and 76 > then there is an intersection of the real parts along the line la and an intersection of the 
imaginary parts along the line lb (in case 7a > and 76 < the lines la and lb are interchanged). 
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see Figure EId. If 7a < and 7b < then the real parts intersect along the both lines la and 
forming a "cluster of shells", while there is no intersections for the imaginary parts, see FigureEt- 
Finally, if 7^ > and 7^ > then there is no intersections for the real parts, while the imaginary 
parts intersect along the both lines la and If,, see Figure Eji. 



4 Example 

As a physical example, we consider propagation of light in a homogeneous non-magnetic crystal 
in the general case when the crystal possesses natural optical activity (chirality) and dichroism 
(absorption) in addition to biaxial birefringence, see Berry and Dennis (2003)| for the general 



formulation. The optical properties of the crystal are characterized by the inverse dielectric tensor 
r]. The vectors of electric field E and displacement D are related as Landau et al. (1984)| 



E = r]B. (48) 

The tensor r/ is described by a non-Hermitian complex matrix. The electric field E and magnetic 
field H in the crystal are determined by Maxwell's equations [Landau et al. (1984) 



--4f • • (-) 

where t is time and c is the speed of light in vacuum. 

A monochromatic plane wave of frequency u that propagates in a direction specified by a real 
unit vector s — (si, S2, S3) has the form 



D(r, t) = D(s) exp iuj \^-^s^ r - t j , H(r, t) = H(s) exp iuj [j^^ r~tj, (50) 

where n(s) is a refractive index, and r = {xi, X3) is the real vector of spatial coordinates. 
Substituting the wave H50|l into Maxwell's equations (|49|) . we find 

H = n[s,r7D], D = -n[s,H], (51) 

where square brackets indicate cross product of vectors [Landau et al. (1984)|. With the vector H 
determined by the first equation of (|51|l , the second equation of l|5Hl yields [Berry and Dennis (2003) 



- [s, [s, r,D(s)]] = r,D(s) ~ s(s^r,D(s)) - ^D(s). (52) 

n^(s) 

Multiplying equation H52|l by the vector s-^ from the left we find that for plane waves the vector 
D is always orthogonal to the direction s, i.e., s-^D(s) — 0. 

Since the quantity S"^J7D(s) is a scalar, we can write (|52|) in the form of an eigenvalue problem 
for the complex non-Hermitian matrix A(s) dependent on the vector of parameters s = (si, S2, S3): 

Au = Au, A(s) = (I-ss^)r7(s), (53) 

where A = n^^, u = D, and I is the identity matrix. Multiplying the matrix A by the vector s 
from the left we conclude that s-^A = 0, i.e., the vector s is the left eigenvector with the eigenvalue 
A = 0. Zero eigenvalue always exists, because det(I — ss^) = 0, if |ls|| 1. 

The matrix A(s) defined by equation H53|) is a product of the matrix I — ss^ and the inverse 
dielectric tensor 77(5). The symmetric part of rj constitutes the anisotropy tensor describing the 
birefringence of the crystal. It is represented by the complex symmetric matrix U, which is 
independent of the vector of parameters s. The antisymmetric part of rj is determined by the 
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Figure 6: Eigenvalue surfaces near a point of weak coupling. 
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optical activity vector g(s) = (31,52,33), describing the chirality (optical activity) of the crystal. 
It is represented by the skew-symmetric matrix 



G = i\ 33 -31 . (54) 
\ -.92 5i / 

The vector g is given by the expression g(s) = 7s, where 7 is the optical activity tensor represented 
by a symmetric complex matrix. Thus, the matrix G(s) depends linearly on the parameters si, 

32, S3. 

In the present formulation, the problem was studied analytically and numerically in [Berry and Dennis (2003) 



Below we present two specific numerical examples in case of a non-diagonal matrix 7, for which the 
structure of singularities was not fully investigated. Unlike Berry and Dennis (2003)| , where the 



reduction to two dimensions was carried out, we work with the three-dimensional form of problem 
(|53|l . Our intention here is to give guidelines for using our theory by means of the relatively simple 
3x3 matrix family, keeping in mind that the main area of applications would be higher dimensional 
problems. 

As a first example, we choose the inverse dielectric tensor in the form 

(55) 
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where S3 ~ -\/l — — s|. The crystal defined by H55I) is dichroic and optically active with the 
non-diagonal matrix 7. When si = and S2 = the spectrum of the matrix A consists of the 
double eigenvalue Aq = 2 and the simple zero eigenvalue. The double eigenvalue possesses the 
eigenvector Uq and associated vector Ui : 

i \ / \ 

(56) 

The eigenvector Vq and associated vector vi corresponding to the double eigenvalue Ao = 2 of the 
adjoint matrix A* are 

vo = I 1 I . VI = I ) . (57) 

\ 1 + 1/2 J \ 1/2-1/4 J 

The vectors Uq, Ui and vq, vi satisfy the normalization and orthogonality conditions Calcu- 
lating the derivatives of the matrix A(si, S2) at the point Sq — (0, 0, 1) we obtain 



— - I i 1 — 




(58) 



Substitution of the derivatives H58|l together with the vectors given by equations l|56(l and (|57|l 
into the formulae H12|l and (|13|l yields the vectors f , g and h, r as 

f=(0,4), g-(-4,0), h=(0,0), r = (-4,0). (59) 

With the vectors H59|) we find from (|21|l and l|22|) the approximations of the eigensurfaces ReA(si, S2) 
and ImA(si, S2) in the vicinity of the point Sq — (0, 0, 1): 



Rc\±=2±\ 2s2 + 2 J sj + sj, lmX±=-2si±\ -2s2 + 2Jsl + sl. (60) 
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Figure 7: Eigensurfaces of the crystal H55(l and their approximations. 

Calculation of the exact solution of the characteristic equation for the matrix A with the inverse 
dielectric tensor rj defined by equation (|^ shows a good agreement of the approximations (|Hn)l 
with the numerical solution, see Figure[7| One can see that the both surfaces of real and imaginary 
parts have a Whitney umbrella singularity at the coupling point; the surfaces self- intersect along 
different rays, which together constitute a straight line when projected on parameter plane. 
As a second numerical example, let us consider the inverse dielectric tensor as 

/ 1 1 \ / 5 4 \ / -si~is2 is3 \ 

t;= 1 +i 5 2 si+is2 -sg . (61) 

V 1 4 y V 4 2 y V -**3 S3 y 

At s = (0, 0, 1), the matrix A has the double eigenvalue Aq = 1 + 5i with two eigenvectors and the 
simple zero eigenvalue. The eigenvectors Ui, U2 of Aq and the eigenvectors Vi, V2 of Aq for the 
adjoint matrix A* are 




These vectors satisfy normalization conditions . Taking derivatives of the matrix A with respect 
to parameters si and S2, where S3 = yl--sf — sf, and using formula . we obtain 

dn = (-2 - 8i, 0), di2 = {6i, -9 - 4i), dai = {-lOi, 7 - 4z), 622 = (0, (63) 

Using (|63|l in formulae H34|l - H36(l . we find approximations for real and imaginary parts of two 
nonzero eigenvalues X± near the point s = (0, 0, 1) as 

Re A± = 1 - si ± V(|c| +Rec)/2, Im A± 5 - 4si - 2s2 ± y^{\c\ -Rec)/2, (64) 

where c = (45 + 8i)sl + 128isiS2 + (-83 + 8i)s|. 

Approximations of eigenvalue surfaces (|64|l and the exact solutions are presented in Figure |S| 
The eigenvalue surfaces have intersections both in (si,S2,ReA) and (si,S2,ImA) spaces. These 
intersections are represented by two different lines la and lb in parameter space, see Figure Eh)- 

5 Conclusion 

A general theory of coupling of eigenvalues of complex matrices smoothly depending on multiple 
real parameters has been presented. Diabolic and exceptional points have been mathematically 
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Figure 8: Eigcnsurfaces of the crystal Htil|) and their approximations. 

described and general formulae for coupling of eigenvalues at these points have been derived. This 
theory gives a clear and complete picture of crossing and avoided crossing of eigenvalues with a 
change of parameters. It has a very broad field of applications since any physical system contains 
parameters. It is important that the presented theory of coupling gives not only qualitative, but 
also quantitative results on eigenvalue surfaces based only on the information at the diabolic and 
exceptional points. This information includes eigenvalues, eigenvectors and associated vectors with 
derivatives of the system matrix taken at the singular points. We emphasize that the developed 
methods provide a firm basis for analysis of spectrum singularities of matrix operators. 
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